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ABSTRACT 


A Navier-Stokes code, FDNS, is used to analyze the complicated internal flowfield of the 
SRM (solid rocket motor) to explore the impacts due to the effects of chemical reaction, particle 
dynamics, and slag accumulation on the solid rocket motor (SRM). The particulate multi-phase 
flowfield with chemical reaction, particle evaporation, combustion, breakup, and agglomeration 
models are included in present study to obtain a better understanding of the SRM design. 

Finite rate chemistry model is applied to simulate the chemical reaction effects. Hermsen 
correlation model is used for the combustion simulation. The evaporation model introduced by 
Spalding is utilized to include the heat transfer from the particulate phase to the gas phase due to 
the evaporation of the particles. A correlation of the minimum particle size for breakup expressed 
in terms of the AI/AI2O3 surface tension and shear force was employed to simulate the breakup of 
particles. It is assumed that the breakup occurs when the Weber number exceeds 6. A simple 
agglomeration model is used to investigate the particle agglomeration. However, due to the large 
computer memory requirement for the agglomeration model, only 2D cases are tested with the 
agglomeration model. The VOF (Volume of Fluid) method is employed to simulate the slag 
buildup in the aft-end cavity of the RSRM. Monte Carlo method is employed to calculate the 
turbulent dispersion effect of the particles. 

The flowfield analysis obtained using the FDNS code in the present research with finite rate 
chemical reaction, particle evaporation, combustion, breakup, agglomeration, and VOF models will 
provide a design guide for the potential improvement of the SRM including the use of materials and 
the shape of nozzle geometry such that a better performance of the SRM can be achieved. The 
simulation of the slag buildup in the aft-end cavity can assist the designer to improve the design of 
the RSRM geometry. 
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1.0 INTRODUCTION 


1.1 THE NATURE OF THE PROBLEM 

It is known that the flowfield of the solid rocket motor is very complicated due to the 
chemical reaction, particle evaporation, combustion, and breakup, and other complex 
characteristics like agglomeration and coalescence etc. Because the distribution of the particles 
affects the performance of the motor, the prediction of the particle effects plays an important role 
for the SRM design. 

Traditionally, metal powders are used in solid propellants for the purpose of increasing the 
motor specific impulse. High density and high heat of reaction are two factors which contribute to 
high impulse. It is known that the simulations of internal flowfields of SRM’s with Al-based 
propellants require complex multi-phase, turbulent, and chemical reacting flow models. On the 
other hand, because of relative velocity and temperature differences between particles and the gas 
flow, inter-phase drag forces and heat transfer exist. Therefore, the effect of particles on the 
flowfield has a significant impact. The evaporation of the A 1 and AI2O3 transfers the mass from 
particulate phase to the gas phase. The combustion of aluminum produces aluminum oxide and 
releases heat and mass. The breakup of the particles affects the flowfield. It is known that a 
recirculation zone near the entry of the aft-dome cavity disturbs the flowfield and increases the 
complexity of interaction between the particle phase and the nozzle inlet section. This interaction 
will determine the slag agglomeration rate which affects the nozzle erosion and motor performance. 
All these impacts should be investigated so that a better performance of the SRM can be achieved. 

To provide design guides for maintaining high performance of the SRM, an accurate 
simulation of the gas-particle interaction is very important. Because of the complex flowfield inside 
the SRM, limited experimental data is available for design purpose. The internal flowfield analysis 
using a CFD (Computational Fluid Dynamics) method can be utilized to obtain a better 
investigation for SRM’s due to the recent progress in computing power. There has been some 
research conducted in the past for the SRM internal flowfield analysis using the CFD method. 
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Madabhushi et al. (Ref. 1) calculated the two-phase aft-dome flowfield of the solid rocket 
motor. The 19-sec bumback configuration was used for the analysis and no particle trajectories in 
the aft-dome cavity was provided. This may not reveal the realistic particle effects for the first 19 
seconds. Since the particles will accumulate on the wall and change the shape of the flow passage. 
Due to the large impact of turbulent particles, a full configuration should be used to include the 
effect of upstream particles. Carrier et al (Ref. 2) investigated the alUminum-oxide-particle field 
within a long-bore SRM with a simple modeling. The Lagrangian particle tracking method was 
used in Ref. 2. Lupoglazoff and Vuillpt (Ref 3) simulated the stability of a 2-D SRM numerically 
by means of the finite volume explicit predictor-corrector McCormack scheme that solves the Euler 
equation. Sabnis et al (Ref. 4) used an Eulerian-Lagrangian two-phase approach to model the 
multi-phase reacting internal flow of a SRM with a metalized propellant. Some other studies in the 
solid propellant rocket motors have been conducted (Refs. 5-8). However, the realistic 
applications including the multi-phase flow with chemical reactions, evaporation, combustion, 
breakup, and agglomeration models were not simulated at the same time to investigate the influence 
of the particles on the SRM. No details about the effects of the recirculation zone on the motor 
performance have been revealed. Due to the complex flowfield at the entry of the cavity, more 
investigation for the aft-dome cavity should be conducted using the numerical approach. Since it is 
very difficult to measure the data for the internal flow of the solid rocket motor. It is found that the 
slag buildup in the cavity causes the oscillated pressure and may affect the motor performance. Not 
many studies have been made to investigate the effects of the slag accumulation. 

To explore the impacts due to the effects of chemical reaction, recirculation zone, slag 
buildup, and particle dynamics, a Navier-Stokes code, FDNS (Refs. 9-11), is used to analyze the 
complicated flowfield for the SRM with chemical reaction, and particle evaporation, combustion, 
breakup, and agglomeration models. 

1.2 OBJECTIVE 

The objective of this project is to develop an advanced particulate multi-phase flow model 
which can simulate the effects of particle dynamics, chemical reaction and hot gas flow turbulence. 
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The inclusion of particle evaporation, combustion, breakup, and agglomeration, particle/gas 
reaction and mass transfer, and slag buildup in modeling the particle dynamics will allow the 
proposed models to realistically simulate the internal flowfield of a solid rocket motor. 
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2.0 CFD METHODOLOGY 


2.1 GOVERNING EQUATIONS 

The FDNS flow solver (Refs. 9-11) provides steady-state and unsteady flowfield solutions by 
solving the following transport equations. The general form of the mass conservation equation, 
Reynolds-averaged Navier-Stokes equations, energy equation and other scalar transport equations 
can be written as: 


dpU 

dt 


d Xi 


pUiU 


+ 


ft - 


dU 

d Xi 


Su (1) 


where p and U = (1, u, v, w, h, k, s and oq) stand for the fluid density and the flow primitive variables 
for the continuity, momentum, energy, turbulence model and species mass-fraction equations 
respectively. This general form of the transport equation has one exception that fluid temperature 
instead of enthalpy is used for the diffusion (heat conduction) terms of the energy equation. The 
source terms SU for the momentum, energy, turbulence model and species mass-fraction equations in 
3-dimensional space x{ can be written for fully conservative form as: 
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where p e = (pj + PtV a represents the effective viscosity which is a sum of the laminar viscosity and 
the turbulence eddy viscosity then divided by the Schmidt or Prandtl number, a. <b, Qt and on are 
the energy dissipation function, heat source and the species source term respectively. Di represents 
the drag forces. Mp denotes the rate of mass addition per unit volume due to inter-phase mass 
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exchange. Vi are the particle velocity components, hv stands for the particle vapor enthalpy. Ur is 
the gas velocity relative to the particle velocity. Pr stands for the turbulence kinetic energy 
production rate which is written as: 


P. = 




d u j 
d x t 


+ 


d x 


\2 


( 0 > 

du k 

i 

\dx k ) 
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The turbulence modeling constants Cj and C 2 are given as 1.43 and 1.92 respectively in the 
standard k-s turbulence model. For the extended model, C 2 = 1.90 is a modeling constant and Cj 
takes a functional form as: 

Ci = 1.15 + 0.25 — 
s 

Turbulence Schmidt numbers for the k- and 8-equation are 1.0 and 1.3 respectively for the 
standard model. For the extended model these two constants are modified to be 0.8927 and 1.15 
respectively. Turbulence Schmidt number for the species mass fraction equation is assumed to be 
0.9. The same value is assumed for the energy equation turbulence Prandtl number (0.7 is used for 
laminar flows). To account for compressibility effect on the turbulence models, two methods of 
model correction are available in the FDNS code. They are: (1) k-equation correction, in which 
the dissipation term in the k-equation is modified by a turbulence Mach number; and (2) s-equation 
correction, in which the in the e-equation is modified by the flow Mach number. An equation of 
state of the following form is used to calculate fluid density and provide closure to the above 
governing equations. 

P 

P RT/ M„ 
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where R, T and Mw stand for the universal gas constant, fluid temperature and the mixture molecular 
weight. The fluid temperature is calculated based on the solution of the fluid enthalpy and the 
JANNAF standard thermodynamics data using a Newton's iteration method for finding the roots of 
the polynomials. 


Particulate-Phase Equations 

The equations that constitute the particle trajectory and temperature history can be written as: 


DVl 

Dt 

Dh, 

Dt 


= [U, -V)lt d 
= Cp c (T aw -T p )lt H 


foeJT* / {pd) p 


where 

Ui = gas velocity 

Vi = particle velocity 

td = particle dynamic relaxation time 

= 4 p p d p l(3Cd p c \u t - 

hp = particle enthalpy 
Tp = particle temperature 
Taw = gas recovery temperature 
tH = particle thermal-equilibrium time 

= (p rf ) ? / [ i2Ar "/ </ ( p ''^)] 

a = Stefan-Boltzmann constant 
= 4.76E-13 BTU/FT 2 -S-°R 
s = particle emissivity = 0.20 ~ 0.31 
f - radiation interchange factor 

Cd and Nu stand for drag coefficient and Nusselt number for heat transfer which are functions 
of Reynolds number and relative Mach number. Typical correlation are described in references 12 
and 13. 
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Turbulent dispersion 

To simulate the turbulent dispersion of particles, the Gaussian probability distribution with 

standard deviation equal to (2k/3)^^ is used for the random turbulent velocity components. The 

turbulent velocity components are thus computed using 
w' = (4k / 2) xl2 erf~ x (lx-]) 

where x is a random variable with uniform probability distribution between 0 and 1 . The generated 
turbulent velocity components are added to the mean velocity field of the continuous phase in 
evaluating the interphase drag force. 

2.2 NUMERICAL ALGORITHM 

The FDNS flow solver is a finite difference method for solving non-linear governing equations 
using non-staggered curvilinear grid system. This code provides multi-zone multi-block options (Ref 
14) for multiple species and finite rate chemistry reacting flow by solving the Navier-Stokes 
equations for the simulation of complex geometry flow problems. A Lagrangian-Eulerian particle 
tracking method is employed in the FDNS to provide effects of momentum and energy exchanges 
between the gas phase and the particle phase. The particle trajectories are calculated using a one-step 
implicit method for several groups of particle sizes by which the drag forces and heat fluxes are then 
coupled with the gas phase equations. A second-order upwind scheme is employed to approximate 
the convection terms. Viscous fluxes and source terms are discritized using second-order central 
difference approximation. The time domain discritization of the present method allows the finite 
difference equations to be arranged into delta form for time-marching integration. Time-centered or 
Euler implicit time-marching schemes are employed for time accurate or steady-state computations 
respectively. A CFL number conditioned non-uniform time marching option can also be used for 
efficient steady-state solutions. The final linearized algebraic equations are solved by iterative point 
relaxation, ADI or L-U matrix solver. The time-marching scheme is described below. For 
convenience, transformed equations (from Xi to £i system with J as the Jacobian of coordinate 
transformation) of Eq. 1 is written as: 
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where F represents convection and diffusion fluxes. First, Eq. 2 is discretized in time with a second- 

order time-centered scheme. That is 

-^{(pr/)"* 1 -{puf} =^(V +1 + V) 

where superscripts n and n+1 represent old and new time levels respectively. If a sub-iteration 
procedure within a time step is applied, the following linearization can be incorporated. 


(pU) t+> ={pU) k + p k AU k 

AC/*+V 

where the superscript k denotes the k-th sub-iteration. With the above approximations, the final form 
of the time-marching scheme can be written as: 

\y (X Y _ ( pu) k - ( pu)" V + s,, m 

JAt \cU J j JAt 2 


The solutions at time level n+1 is then updated by: 

U n+] = U k+l =U k +AU k 


When k = 1 is selected, a non-iterative time-marching scheme with a multi-corrector solution method 
can provide time accurate solutions for unsteady flow problems. The pressure based multi-corrector 
solution method is formulated using simplified perturbed momentum and continuity equations. The 
simplified velocity correction equation can be written as: 


d pu t 
dt 


« -VP 


or, in discrete form, 

At , 

u^-p—VP' and P n+l =P n +P' (3) 
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where J3 represents a pressure relaxation parameter (typically 10). The velocity and density fields in 
the continuity equation are then perturbed to form a correction equation. Higher order terms are 
neglected. That is, 

(4) 

Substituting Eq. 3 into 4 and letting p' = P'/RT, the following all speed pressure correction equation 
is obtained. 

+ = -V{puX (5) 

RT dt VRT J ’ \dt) yH 

To provide smooth shock solutions the adaptive dissipation terms based on the pressure field is 
added to the right hand side of Eq. 5. Once solution of Eq. 5 is obtained, the velocity and pressure 
fields are updated using Eq. 3. The density field is then updated through the equation of state. The 
temperature field can also be modified by using a perturbed temperature correction equation. The 
entire corrector step is repeated 3 or 4 times such that the mass conservation condition is enforced 
before marching to the next time level. 
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3.0 PROPOSED MODELS 
3.1 FINITE RATE CHEMISTRY MODEL 

For gas-phase chemical reaction modeling, a general system of chemical reactions can be 
written in terms of its stoichiometric coefficients (vy and vy' ) and the i-th chemical species name 
(Mj) of the j-th reaction as 

X VtjMi = ^ Vij'Mi 

i i 

The net rate of change in the molar concentration of species i due to reactions j, Xy, can be 
written as: 

X l/ = (v li '-v l ,)(Kj e n(pa l /M m ) V,i 

and the species production rate, ©i, (in terms of mass fraction) is calculated by summing over all 

reactions. 

(Dj — h/f TX,, 
j 

where 


M wi = 

molecular weight of species i 

(Xj = 

mass fraction of species i 

P 

fluid density 

K 5 = 

forward rate of reaction j 

K b j = 

backward rate of reaction j = Ky/K e j 

Kei = 

equilibrium constant 

= 

(MRTf^EXP^f: v, -f t v,)} 

fi = 

Gibbs free energy of species i 
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K f = A T b EXP(-E/RT) 

Finally, the species continuity equations can be written as: 


where c a (assumed to be 0.9) represents the schmidt number for turbulent diffusion. A penalty 
function is employed to ensure the basic element conservation constraints at the end of every time 
marching step. This is a crucial requirement for the numerical stability and accuracy of a CFD 
combustion model. It is accomplished by limiting the allowable changes in species concentrations, 
which are the solutions of the species continuity equations, for each time step such that the species 
mass fractions are well bounded within physical limits. The resulting limited changes are adjusted 
so that they are proportional to the species source terms. 

3.2 COMBUSTION MODEL 

The combustion model using Hermsen correlation (Ref. 15) is employed. The complete 


description of this correlation is described in the following: 


d( m 4i) x 

dt 2 Pa ' 



n 


where the exponent n is 1.8, Dp is the particle diameter, pAl is the density of aluminum, and the 
burning rate constant k (cm 1 8 /s) in the above equation is computed from 
K = 8.3314 x 10- 5 Ak°- 9 Pc 0 - 27 Sh/2 

where Pc is the chamber pressure (psia), Sh is the Sherwood number based on particle diameter, 
and Ak is a measure of the availability of oxidizing species and is computed using 
A k = 1 00Z Xf ; /' = COi, HiO, O2, OH, O 

i 
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The Sherwood number was computed using empirical correlation relating it to the Reynolds 
number, 

Sh = 2.0 + 0.6 Re D ° 5 Pr 0687 ; Re < 278.92 

p p 5 p 

Sh p = 2.0 + 0.6 Re/ 5 Pr° 687 ;K& p > 278.92 

The calculated mass loss (nui) of particles due to the combustion at each time step is used to 
compute the new particle size. The combustion model is used for the A1 particles only. A1 particle 
starts to bum when its surface temperature exceeds its melting point. The diameter of A1 particle 
reduces when it bums and the aluminum oxide is formed on the particle surface. It is assumed that 
the combustion stops when the particle diameter is reduced to 50pm. The aluminum becomes 
aluminum oxide with diameter of 50pm when the combustion stops. 

The calculated mass loss at each time step is used to compute the new particle size. For this 
study, D 43 , the mean diameter of particle at the nozzle exit is used to determine when the 
combustion is completed. 

D 43 = 3.6304 Dt 0 - 2932 (1 - e - 0000816sPx ) 
where 

D 43 = mass-weighted average diameter, pm. 

Dt = nozzle throat diameter, in. 

8 = Aluminum concentration in chamber, g-mole/lOOg. 

P = chamber pressure, psia. 
x = average chamber residence time, msec. 

3.3 EVAPORATION MODEL 

For gas temperature exceeding the metal particle boiling point (e.g. 2750 °K for aluminum 
oxide particle under 1 atm pressure), inter-phase mass transfer can occur through evaporation or 
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dissociation processes. A commonly adopted treatment for evaporating particles can be 
summarized as follows: 

The mass transfer rate from the discrete particle to the primary gas phase associated with the 
evaporation is given by Ref. 16: 

Mp = (7i/2) p P d P 2 (Ddp/Dt) = 27td P (K/Cp) ln(l +B) 
where 

Mp = particle mass transfer rate due to evaporation 
K = thermal conductivity of continuous phase 
Cp = specific heat of continuous phase 
B = mass transfer number, Cp,v(T« - Tw)/L 
Cp,v - specific heat of vapor 
Too = ambient temperature of continuous phase 
Tw = particle surface temperature 
L = latent heat of vaporization at temperature Tw 

The evaporation rate Mp will be reduced when the aluminum oxide shell forms on the particle 
surface due to the combustion. A correction factor depending on the thickness of the shell is 
applied to include the effect of the oxide shell on the evaporation. The evaporation model is used 
for AI2O3 particles due to the lack of latent heat of vaporization of Al. 

3.4 BREAKUP MODEL 

This model uses Weber number to determine if the breakup of the particle occurs. After the 
particle temperature reaches the melting point of AI2O3 and AL, the liquid AI2O3 and AL could 
break up based on the Weber number (Ref. 17): 

We = Dp pg (Ug - Up) 2 / a 
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where 


a = the surface tension, 

We = the critical Weber number for particle breakup. 

Dp = diameter of particle/agglomerate 
pg = density of gas phase 
Ug = velocity of gas phase 
Up = velocity of particle/agglomerate 

The surface tension of AI2O3 at 2,300°K is 0.69 N/m (summarized in Ref. 18). The surface 
tension of AL is between 0.85 and 0.90 N/m in the temperature range of 970 to 1,020 °K 
(tabulated in Ref. 18 and 19). For this study, 0.75N/m for the temperature less than 2,300°K and 
0.65 N/m for the temperature less than 2,300°K are assumed to determine the surface tension for 
AI2O3 particles/agglomerates. Linear interpolation/extrapolation is used to compute the surface 
tension of AL particles based on the known data. It is known that for low Weber numbers, 
droplets are spherical. They tend to distort when the Weber number exceeds about 4. The 
distortion increases with increasing Weber number until breakup occurs in the range of We = 20 ~ 
30. We = 6 is used for current study. 

The breakup time of a droplet suddenly exposed to a gas stream has been approximated as 
(Ref. 20): 


d 


ag 


2(u g -u ag ) 




1/2 


ty, is approximately 0.1 ms, which is very consistent with the breakup events observed in the 
experiments conducted in Ref. 20. 


3.5 AGGLOMERATION MODEL 
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The transient simulation is needed for the particle agglomeration. Due to the large computer 
memory requirement for using the Lagramgian method to simulate the particle agglomeration, 
assumptions for the agglomeration model are needed: 

(1) Agglomeration model is turned on after the statistic particle trajectories are determined using 
the combustion, evaporation, and breakup models. 

(2) Agglomeration occurs only when Weber number is less than 5 and the distance between two 
particle centers are shorter than 3 quarters of the sum of their radii. 

(3) The particles merged together and the new diameter is computed based on the total mass of the 
merged particles. 

Summary of the models 

Figure 1 illustrates the models of combustion, evaporation, breakup, and agglomeration used 
in present study. The aluminum particles (Al) with 150pm combust and forms the aluminum oxide 
(AI2O3) while the diameter is reduced due to the combustion and evaporation. The combustion 
model is turned off when the diameter of the particle is less than 50pm. The particle is considered as 
aluminum oxide when its diameter is less than 50pm. The size of the aluminum oxide keeps 
shrinking because of evaporation and breakup. When the diameter of the particle reduces to 5pm, 
which is assumed as the maximum diameter of the smoke, the size will not change any more. The 
particle agglomeration occurs when the conditions described in 3.5 are satisfied. 

3.6 VOF MODEL 

The VOF (Volume Of Fluid, Refs. 21 and 22) method is used to account for the effect of slag 
buildup in the aft end region. VOF method can be used to predict the sloshing dynamics, in 
response to the flight dynamics and local acceleration of the slag in the aft end region. 
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Traditionally, VOF methods are mainly developed and used for low-speed flows such that 
incompressibility can be assumed. The incompressible flow assumption has limited their capability. 
To generalize, the present formulation is based on compressible flow governing equations. The 
forms of the equations are then continuously reduced to their incompressible forms according to the 
local flow conditions and the VOF solutions. To illustrate this, a general transport equation can be 
written as: 


dps 

a + dx t 



And, the VOF transport equation is given below: 


da t \aa 


da 


( 6 ) 


where a = 1 stands for liquid and a = 0 is for gas. The interface is located atl>a>0. For a 
given solution of a field, equation (6) can be recast as: 



a > 0.01 
a < 0.01 


for compressible gas 
for incompressible gas 


and 

p m =Max[p s ,ap l ) 

where p g and pi denote gas and liquid density respectively, ug represents the grid speed 
components used to simulate moving domain effects. The numerical accuracy of the VOF method 
depends highly on the interface resolution. To prevent the solution from becoming too smearing 
due to numerical diffusion, a compression procedure is developed to perform VOF interface 
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rescaling such that the total volume within the interface (0.1 <a<0.9)is kept constant through 
out the computation. The interface a solution compression procedure is expressed as: 


= Max\o , M'k[ 1 , 0.5 + f(a M - 0.5)]} 
and 


/= 


( Interface volume) 

(interface volume ) 

The surface tension forces in the continuum surface force model is formulated as continuous 


initial 


body forces across the interface. These forces can be written as: 

F, =-a (vnJa Jt 


F,= 


- CF^ V Tl)jCC ) 


+ 


a. 


K y J 


for 2D axisymmetric only 

F z =~ a^V n^a z , for 3D case only 

where 

= surface tension constant 


A A A A 

Vn = a xx 4 -ayy +a Z2 


a is 0.5 for the free surface. The VOF method is used to represent the tracking of the free 
surface between the liquid and gas phase. 
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4.0 NUMERICAL RESULTS 

To test the proposed models, a simple 2-D motor, the 2D ASRM, the 3D ASRM, and a 2D 
RSRM are used for present study. The computational results show that a recirculation zone exists 
at the entry of the aft-dome cavity. The particle impingement could cause erosion and damage the 
nozzle wall. The particles may accumulate in the impingement area and change the wall shape and 
affect the performance of the motor. The accumulation of the slag in the aft-end cavity may affect 
the performance of the solid rocket motor. The flowfield is disturbed by the particles and the slag. 
The pressure jump due to the slag accumulation is captured in the present study. The obtained 
results are comparable to the known data. The phenomenon is reasonable based on the physical 
point of view and available experimental data. 

4.1 2-D NOZZLE 

A simple 2-D axisymmetric solid rocket motor configuration (Fig. 2) of 8 lx 31 grid points 
generated with algebraic method is used for present study. Launch conditions with lg gravity force 
at 1 atm are considered. The gas phase flow with chemical reaction using 7 species (H20, 02, H2, 
O, H, OH, and N2 ) with 9-step chemical reaction is simulated for the 2-D SRM geometry with and 
without three groups of particles respectively. 

(1) Without particles 

The Mach number distributions without chemical reaction and particle effects are shown in 
Fig. 3. The flow accelerates because of the mass flow issuing from the propellant surface. It is 
seen that Mach number reaches 1.0 at the nozzle throat. The maximum Mach number of 3.91 at 
the nozzle exit is reached. Figure 4 shows the Mach number contributions with the simulation of 
chemical reaction. The maximum Mach number is reduced to 3.642 due to the effect of the heat. 
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The flowfield is smooth without the effects of particles. Figure 5 illustrates the pressure contours 
without the effects of chemical reaction and particles. The constant pressure contours in the 
chamber at each cross station is obtained. The maximum chamber pressure is 929.6psi. The 
pressure contours with the effect of chemical reaction are shown in Fig. 6. The maximum chamber 
pressure is increased to 943.3psi due to the effect of chemical reaction. The obtained pressure 
distributions are also comparable to that of the regular SRM. Lower temperature is obtained 
(shown in Fig. 7) for the case without chemical reaction. While higher temperature due to the 
effect of the chemical reaction is shown in Fig. 8. This phenomenon is reasonable. 

(2) With particles 

Based on the observation by Price (Ref. 23), three major particle groups injecting from the 
grain surface are assumed to simulate the particulate phase. Particles with diameter of 150 or 
175pm represent the A1 particles which bums on the grain surface. The second group of particles 
with diameter of 50pm is assumed to represent the aluminum oxide (AI2O3) which is the product 
due to the combustion of A1 on the propellant surface. The third group of particles is the “smoke” 
which has maximum diameter of 5pm. It is assumed that the “smoke” does not evaporate, 
combust, break up, and agglomerate. All three groups are injected from the propellant grain 
surface vertically into the chamber. In order to make comparison, a test case without 
combustion/evaporation/breakup/agglomeration models is also conducted. 

Figure 9 illustrates the trajectories of three groups of particles. It is seen that the size of each 
particle group is constant. The models of combustion/evaporation/breakup are then used for the 
test. Figure 10 illustrates the particle trajectories due to the effects of combustion, evaporation, 
and breakup models. The aluminum particles (Al) with diameter of 150pm combust and form the 
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aluminum oxide (AI2O3) while the diameter is reduced due to the combustion and evaporation. 
Because of the high temperature in the chamber, A 1 particle becomes AI2O3 within a short period of 
time after leaving the propellant surface. The AI2O3 oxide with diameter larger than 5pm does not 
bum but will evaporate and may break up. The burning time of the A 1 particle is longer for those 
particles coming out of the propellant surface in lower-temperature areas near the nozzle. It is 
noted that the particle trajectories are not plotted at each time step. The AI2O3 particle is 
accelerated when it approaches the nozzle. The relative speed of the particle and the flow and the 
change of particle surface tension due to the temperature difference allow particles to break. It is 
seen that most of the AI2O3 particles break up at the nozzle throat. This phenomenon is identical to 
the experimental data shown in Ref. 17. The breakup model requires large agglomerates formed 
before the agglomerates enter the nozzle. Experimental data of Ref. 17 demonstrates that the 
breakup occurs mostly near the throat. The computed results shown in Fig. 1 1 using current 
models match the experimental conclusion in Ref. 17. The average of the particle diameter reduced 
to less than 10pm after the nozzle throat. The agglomeration model is tested and shown in Fig, 12. 
The AI2O3 particle agglomerates are based on the Weber number and the distance between the 
particles. The particles merge together and become larger when they agglomerate. Figure 12 
demonstrates the computational results. It is noted that only few groups are shown in Fig. 12, 
because it is very difficult to show the particle trajectories for particle agglomeration. The mean 
particle size distributions shown in Fig. 13 are used to estimate the averaged particle size. It is seen 
that the size of the particles becomes less than 10pm after they pass the nozzle throat. Figure 14 
shows the Mach number contours due to the effects of chemical reaction and particles. The effects 
of particles on the flow field is clearly seen by comparing Fig. 14 and Fig. 3. Mach number is 
reduced due to the drag of particles. Apparently, the flowfield is disturbed by the particles, i.e., the 
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flow speed is reduced due to the effect of particle concentrations away from the center line. The 
maximum Mach number at the nozzle exit reduced to 3.38 from 3.91. Pressure contours are shown 
in Fig. 15. Higher pressure (99.6atm) is obtained with the particle effects. This is one of the 
reasons why A 1 is used for one of the propellant ingredients. Higher pressure in the chamber can 
provide more power for the SRM. Figure 16 shows the temperature contours. The temperature is 
higher due to the heat release from the particles. The increase of the temperature due to the effects 
of particles is about 163°R. The temperature contours also depend on the particle concentrations. 
A constant temperature contour in the chamber exists due to the effects of particles. This constant 
temperature contour line across the chamber should not exist if a more accurate analysis using 
various initial particle size distributions on the propellant surface is performed. 

4.2 2-D ASRM 

A simplified 2-D ASRM geometry with axisymmetric flow assumption at launch condition (1-g 
gravity is assumed) is used for the present study. Figure 17 shows the ASRM configuration. Figure 
18 shows the grid system of the 2-D ASRM geometry including the front grain port, the inhibitors, 
and the aft-end cavity. Algebraic method was used to generate the grids. The geometry and mass 
flow rate information are provided by the NASA Marshal Space Flight Center. The gas phase flow 
with chemical reaction using 12 species (H20, 02, H2, O, H, OH, CO, C02, CL, CL2, HCL, and 
N2) and 18-step chemical reactions is simulated for the 2D ASRM geometry. 

(1) Without particles 

The Mach number distributions without chemical reaction and particles are shown in Fig. 19. 
The maximum Mach number is about 3.406 at the nozzle exit. The effect of chemical reaction on 
the Mach number is shown in fig. 20. The Mach number is reduced to 3.15. The pressure contours 
are shown in Fig. 21. The maximum pressure is 924psi in the combustion chamber. The 
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temperature is about 6,314°K in the combustion chamber. The velocity vectors near the aft-end 
cavity are shown in Fig. 23. The recirculation zone near the entry of the cavity is captured. It is 
believed that the impingement of the particles in this area may cause the erosion, damage the nozzle 
surface, and affect the performance of the motor. 

(2) With particles 

The particle trajectories without combustion/evaporation/breakup/agglomeration models near 
the front grain port, the first inhibitor, and the nozzle are shown in figures 24, 25, and 26 
respectively. It is noted that the particles numbers are less and the particle sizes are smaller in fig. 
26 in order to have a clear picture. It shows that the upstream particles affect the downstream 
flowfield. The flow speed in the cavity is very low. Due to the pressure gradient, a recirculation 
zone is formed near the entry of the cavity. The flow recirculation is the mechanism to change the 
particle trajectories. The particles are turned due to the recirculation and travel along the grain 
surface through the cavity and into the nozzle. Because of the unstable nature of the flow in the 
cavity, the steady solution can not be obtained for the flowfield inside the cavity. These show the 
complicated flow phenomena in the area near the cavity and the nozzle. 

Figure 27 illustrates the mean particle size distributions near the nozzle using combustion, 
evaporation, breakup, and agglomeration models. It is seen that the maximum particle size is about 
48pm. The particles size become smaller due to the breakup at the nozzle throat. The maximum 
particle size at the nozzle is about 20pm. The larger particles come from the grain surface near the 
entry of the nozzle. A more accurate set of initial conditions of the grain surface is required to 
obtain a better simulation. Figure 28 shows the Mach number distributions near the nozzle. 
Apparently, the flowfield is disturbed by the particles, i.e., Mach number distributions in the nozzle 
are deformed and the flow speed is reduced due to the effect of particle concentrations away from 


24 



the center line. The Mach number is less than 0.35 in the chamber and ranges from 2.65 to 3.01 at 
the nozzle exit 

4.3 3DASRM 

The 3D ASRM geometry with axisymmetric flow assumption at launch condition (with 1-g 
gravity) is used for the present study. Figure 29(a) shows the grid system of the 3D ASRM geometry 
including the front grain port, the inhibitors, and the aft-end cavity. Only one eleventh of the whole 
geometry is simulated due to the symmetric geometry. Figure 29(b) illustrates the domain near the 
front grain port and the first inhibitor. Figure 29(c) illustrates the domain near the aft-dome cavity 
and the submerged nozzle. It notes that Figure 29 is a rough schematic illustration showing how the 
grids are distributed. An algebraic method was used to generate the grids. 64,152 grid points are 
used for the whole grid system. The gas phase flow with chemical reaction using 12 species (H20, 
02, H2, O, H, OH, CO, C02, CL, CL2, HCL, and N2) and 18-step chemical reactions is simulated 
for the 3D ASRM geometry with and without particles respectively. 

(1) Without particles 

The test is first conducted for the case with chemical reaction and without particles. The Mach 
number distributions are shown in Figure 30. It is seen that Mach number reaches 1.0 at the nozzle 
throat. Very low Mach number is obtained near the front grain port. The flowfield near the end of 
the front grain port is disturbed due to the flow coming out of the grain slots. The Mach number is 
less than 5.6 in the combustion chamber. The flow accelerates and reaches the maximum Mach 
number of 3.508 at the nozzle exit. The pressure contours are shown in Fig. 31. More complicated 
pressure distributions near the end of the front grain port and the entry of the aft-end cavity are 
captured. The geometry change in these areas causes the flow disturbance. The simulated chamber 
pressure is about l,200psi. Figure 32 shows the temperature distributions. The maximum 
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temperature is about 6,3 14°R in the combustion chamber and reduced to about 2,700°R at the nozzle 
exit. Figure 33 demonstrates the flowfield near the aft-end cavity. A strong recirculation zone is 
captured. The strong impingement will cause the damage to this area and affect the motor 
performance. 

(2) With particles 

The particle trajectories are shown in fig. 34 for the whole domain. Figure 35 illustrates the 
particle trajectories and velocity vectors near the front grain port. The particles are driven by the 
flow field. It is seen that the particles do not accelerate in the low speed zone near the front of the 
grain port. The particles accelerate and come out of the grain port slot as shown in Fig. 35(a). This 
is demonstrated by the velocity vectors shown in Fig. 35(b) and (c). Figure 36(a) illustrates the 
particle trajectories near the first inhibitor. The particles coming out of the inhibitor are affected by 
the flow field in the chamber as shown. Figure 36(b) shows the velocity vectors near the first 
inhibitor. It is seen that the flowfield of the combustion chamber is affected by the flow coming out 
of the inhibitor. 

Figure 37 shows the complicated unsteady flow characteristics in the aft-end cavity. The 
particles may become slag and accumulate in the cavity. This will affect the motor performance. The 
flow speed in the aft-end cavity is very low. Due to the pressure gradient, a recirculation zone is 
formed at the entry of the cavity. The flow recirculation is the mechanism to change the particle 
trajectories. The particles are turned due to the recirculation and travel along the grain surface 
through the cavity and into the nozzle. Some particles go into the cavity after impinging to the 
nozzle wall. This impingement could damage the nozzle wall. Also, an agglomerate could occur on 
the wall and affect the motor performance. Because of the unstable nature of the flow in the cavity, 


26 



study. No evaporation/combustion/breakup/agglomeration models are applied. Only the 
qualitative analysis is attempted. 

(1) Without particles 

Figures 42, 43, and 44 show the contours of Mach number, pressure, and temperature, 
respectively. Smooth distributions are obtained without the effects of particles. Maximum pressure 
is 42.15atm (619.6psi). Maximum temperature is 6,137°R. These are based on the data provided 
by NASA Marshall Space Flight Center. Figure 45 illustrates the velocity vectors in the aft-end 
cavity without particles. The vortex and an impinging stagnation point on the wall are predicted. 
Due to the large flow velocity difference between the chamber and the cavity, a vortex should exist. 
The predictions of the flowfield near the cavity are reasonable based on the physical phenomena. 
The impingement may damage the material of the wall and affect the performance of the motor. 

(2) With particles 

Figure 46 shows that the molten particles may enter the cavity and accumulate on the wall. 
The slag accumulation depends on the temperature, the vortex strength, and the flight angle. The 
shape of the aft-end cavity dominates the shape and strength of the vortex. Apparently, the slag 
accumulation exist in the aft-end cavity as soon as the propellant starts to bum in the aft-end cavity. 
In order to investigate the effects of the particle size on the slag accumulation, different particle 
diameters are used for this purpose. Figure 47 demonstrates that the slag flow rate entering the 
cavity depends on the particle size. It shows that an efficient particle combustion results in less slag 
accumulation in the aft-end cavity. Also, Fig. 47 shows that a two-way coupling to include the 
interaction of gas and particles is needed for a better analysis. The effects of the particles on the 
Mach number and temperature contours are shown in figures 48 and 49, respectively. The 
contours are affected strongly by the particle concentration. The slag accumulation in the aft-end 
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cavity will change the particle concentration and affects the flowfield. Figure 50 illustrates this 
phenomenon. The slag accumulation in the aft-end cavity at different times shown in Fig. 50 clearly 
demonstrates the effects of the slag on the flowfield. Figure 50(1) shows the initial conditions of 
the assumed slag accumulation and the particle trajectories near the cavity at t = 0 second. It 
shows that the particles impinge on the wall near the entry of the lower surface of the cavity. The 
slag should accumulate in this area. Figure 50(2) (at t = 0.000444 sec) demonstrates this predicted 
phenomenon. It is seen that particles stick on the wall and form the slag layer. The flowfield will 
then be disturbed due to this slag accumulation. Figure 50(3) illustrates that more slag accumulated 
in the cavity at t = 0.0598 sec. The slag accumulated on the upper surface of the cavity. It also 
shows that the slag buildup in the cavity apparently changes the flowfield of the aft-end cavity. 
Figure 50(4) shows that the slag starts to flow out of the cavity along the wall at t = 0.3482 sec. 
The particle concentration near the entry of the nozzle is changed by the slag. The slag in the 
cavity moves due to the oscillated acceleration. More slag comes out of the cavity as shown in 
figures 50(5)-(8). The slag grows by the merge of the particles coming from the chamber. It is 
shown that the particle trajectories differ due to the effects of the slag accumulation. The slag on 
the upper surface of the cavity moves due to the effects of the flowfield and the acceleration of the 
flight (also see figures 50(5)-(8)). The slag will finally enter the nozzle and affect the performance 
of the solid rocket motor. Figure 5 1 shows the pressure history. The slag accumulation increases 
the pressure in the cavity . The flight test shows the same trend. It is noted that the predicted 
results are based on a fixed 2D geometry with particles of 100pm. A better comparison can be 
made if more information from the test flight can be employed for a 3D test. 
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5.0 CONCLUSIONS AND DISCUSSION 


The obtained computational results using FDNS with the proposed models demonstrate the 
complex internal flowfield of the solid rocket motor. The prediction of the recirculation zone in 
the aft-dome cavity, the particle impingement on the wall, the effects of the particles on the 
flowfield and the slag buildup in the aft-end cavity are very crucial for the improvement of the solid 
rocket motor performance. The predicted results are comparable to the known design values and 
the flowfield is reasonable based on the physical point of view. The flowfield analysis using the 
FDNS code in the present study can provide a design guidance for the solid rocket motor. The 
obtained results can provide the designers a basic guide line for the use of materials and the design 
of the geometry. The slag accumulation analysis plays an important role for the SRM design, since 
the distribution of the slag changes the flowfield at the entry of the nozzle and affects performance 
of the motor. A better performance of the solid rocket motor can be achieved by modifying the 
geometry of the aft-end cavity using the CFD method to prevent the formation of vortex and slag 
accumulation in the aft-end cavity. The geometry of the propellant grain can also be improved 
using CFD method to increase the combustion efficiency. 
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1 Illustration of combustion, evaporation, breakup, and agglomeration models. 
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Figure 8 Temperature contours of 2D nozzle, with chemical reaction, no particles, v 
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Figure 9 Particle trajectories without combustion/evaporation/breakup models. 
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Figure 10 Particle trajectories with combustion/evaporation/breakup models. 
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Figure 12 Particle trajectories with combustion/evaporation/breakup/agglomeration 
models. 
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Figure 14 Mach number contours of 2D nozzle, 
with chemical reaction and particles. 
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Figure 18 2D Gad system of the ASRM. 
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Figure 22 Temperature contours of 2D ASRM, with chemical reaction, no particles. 
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Figure 27 Mean particle size distributions near the nozzle. 
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Figure 29 3D Grid system of ASRM. 





Figure 30 Mach number contours of 3D ASRM, 
with chemical reaction, no particles. 
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Figure 41 2D Grid system of RSRM configuration. 



Figure 42 Mach number contours of RSRM, no particles. 
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Figure 48 Mach number contours of RSRM* with particles. 
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